ASYMPTOTIC PRESERVING IMPLICIT-EXPLICIT RUNGE-KUTTA 
METHODS FOR NON LINEAR KINETIC EQUATIONS 

GIACOMO DIMARCO* AND LORENZO PARESCfflt 

Abstract. We discuss Implicit-Explicit (IMEX) Runge Kutta methods which are particularly adapted to stiff 
kinetic equations of Boltzmann type. We consider both the case of easy invertible collision operators and the 
challenging case of Boltzmann collision operators. We give sufficient conditions in order that such methods are 
asymptotic preserving and asymptotically accurate. Their monotonicity properties are also studied. In the case of 
the Boltzmann operator, the methods are based on the introduction of a penalization technique for the collision 
integral. This reformulation of the collision operator permits to construct penalized IMEX schemes which work 
uniformly for a wide range of relaxation times avoiding the expensive implicit resolution of the collision operator. 
Finally we show some numerical results which confirm the theoretical analysis. 
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1. Introduction. The numerical solution of Boltzmann-type equations close to fluid regimes 
represents a real challenge for numerical methods. In these regimes, in fact, the intermolecular 
collision rate grows exponentially and the collisional time becomes very small. On the other hand, 
the actual time scale for evolution is the fluid dynamic time scale, which can be much larger 
than the collisional time. A non dimensional measure of the importance of collision is given by 
the Knudsen number which is large in the rarefied regions and small in the fluid ones. Standard 
computational approaches lose their efficiency due to the necessity of using very small time steps 
in deterministic schemes or, equivalently, a large number of collisions in probabilistic approaches. 
Unfortunately the use of implicit solvers originates a prohibitive computational cost due to the 
high dimensionality and the nonlinearity of the collision operator. 

Several authors have tackled the above problem for the Boltzmann equation in the recent past 
(see [2] [IIS HH [H] and the references therein) . In summary, the possible approaches that permit 
to overcome such a difficulty can be subdivided into two main classes. Domain decomposition 
strategies and asymptotic preserving schemes. The first class of methods permits to avoid the 
problem of very small Knudsen number by identifying the regions where it is possible to use the 
reduced fluid model and the regions where the full kinetic model must be solved. The literature is 
this direction has a long history we recall here references [BJ j33] and some recent works by Degond 
and coauthors [T^l [T3] • A closely related research approach combines stochastic and deterministic 
solvers in the different regions by originating hybrid methods |15i 114] . 

Concerning the asymptotic preserving strategies, these techniques permit to solve the full 
problem in the entire domain for all choices of time steps and Knudsen numbers. Along this 
direction we quote the pioneering papers by Coron and Perthame for the BGK model and by 
Gabetta, Pareschi and Toscani |21) where an explicit exponential technique for the full Boltzmann 
equation capable to avoid the costly inversion of the collision term has been developed. More 
recently several improvements to the above approaches has been presented, we refer to pi 1191 [T5] . 

In the present paper we develop Implicit-Explicit Runge-Kutta methods [TJ [8j [30] which are 
particularly efficient for stiff non linear kinetic equations. In particular we generalize the approach 
recently introduced in [191117] . First we consider the case where the implicit inversion of the collision 
term does not represent a problem, like for example the case of simple BGK operators. Asymptotic 
preservation properties and monotonicity are carefully studied and analyzed. Subsequently we deal 
with the challenging case of the full Boltzmann equation. To this aim, following |19] we introduce 
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a penalization strategy based on a decomposition of the gain term of the collision operator into 
an equilibrium and a non equilibrium part. This permits to derive new penalized IMEX schemes 
which keep the good asymptotic preservation properties of standard IMEX schemes by avoiding 
the costly inversion of the collision term. Similarly to [301 116] , we derive sufficient conditions for 
asymptotic preservation and asymptotic accuracy. Monotonicity property are also considered and 
discussed. 

We emphasize that the penalized IMEX schemes here developed apply to any large system of 
stiff ordinary differential equations of the form 

Y' = F(Y) + - £ R{Y), Y(t )=Y , (1.1) 

where e > is a small parameter, Y, F(Y) 6 M. N and the non-linear operator R(Y) is a dissipative 
relaxation operator [10] . Such operator is endowed with a n x TV matrix Q of rank n < N such 
that QR(Y) = 0, V Y. This gives a vector of n conserved quantities y = QY. Solutions which 
belong to the kernel of the operator R(Y) = are uniquely determined by the conserved quantities 
Y = E(y) and characterize the manifold of local equilibria. 

The methods here studied are based on the following decomposition 

R(Y) = N(Y) + L(Y), (1.2) 

where N(Y) represents the non-dissipative non-linear part and L(Y) is a linear term such that 
L(Y) — implies Y — E(y). For example L(Y) = A(E(y) — Y) where A > is an estimate of 
the Jacobian of R evaluated at equilibrium. Note that, at variance with standard linearization 
techniques which operate on the short time scale, the operator is linearized on the asymptotically 
large time scale. This decomposition permits to apply IMEX techniques which are implicit in the 
linear part and explicit in the non-linear part. The use of such techniques, as we will see, permits 
to achieve unconditionally stable and asymptotic preserving methods at the cost of an explicit 
scheme. 

The rest of the paper is organized as follows. First in Section 2 we recall some basic aspects 
on kinetic equations and their fluid dynamic limits. The notion of asymptotic preservation is also 
introduced. In Section 3, we consider standard IMEX Runge-Kutta schemes applied to kinetic 
equations and derive conditions for asymptotic preservation and asymptotic accuracy. Next, in 
section 4 we introduce the penalized IMEX schemes for the full Boltzmann model. We analyze their 
asymptotic preservation and asymptotic accuracy properties. Monotonicity in the homogeneous 
case is also discussed. Finally some numerical examples of schemes up to third order and conclusions 
are reported in Sections 5 and 6. 

2. The Boltzmann equation and related kinetic equations. We consider kinetic equa- 
tions of the form [5] 

d t f + v-V x f = Q(f), (2.1) 

with initial data f{x, v, i)|t=o = Here f(x,v,t) is a non negative function describing the 

time evolution of the distribution of particles with velocity v € R 3 and position x £ C M. dn: 
at time t > 0. For notation simplicity in the sequel we will omit the dependence of / from the 
independent variables x, v, t unless strictly necessary. The operator Q(f) characterizes the particles 
interactions and in the case of the quadratic Boltzmann collision operator of rarefied gas dynamics 
it reads 

QbU)= I B(\v-v*\,n)[f(v')f(v'*)-f(v)f(v*)]dv*dn (2.2) 
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where 

v' = v + ^(v - v*) + hv - v*\n, u» = v + -(v - u*) - i|t> - u*|n, (2.3) 

and -B(|f — i>*|, n) is a nonnegative collision kernel characterizing the details of the collision. It is 
described by the following equation 

B(\v -v*\,n) = a ( ^ — ^4 • n ) \v ~ v*| 7 , 

with 7 G [0,3). The case 7 = 1 is referred to as hard spheres case, whereas the simplified situation 
7 = 0, is referred to as Maxwell case. Note that in most applications the angle dependence is 
ignored and a is assumed constant. 

The operator Q(f) is such that the local conservation properties are satisfied 



<P(v)Q(f)dv=:{4>Q(f)}=0 (2.4) 
where <j>(v) — ( 1, v, ^j-J are the collision invariants. In addition it satisfies the entropy inequality 
:H{f)= I Q(/)log/d«<0, H(f)= [ flogfdv. (2.5) 



d_ 

dt 



The functions such that Q(f) = are the local Maxwellian equilibrium functions 



P - " |2 



M[/] = M(p, u, T) = exp l^-^r^ ) , (2.6) 

where p, u, T are the density, mean velocity and temperature of the gas in the apposition and at 
time t defined as 

(p,pu,E) T = {4>f), T = ±(E - p\u\ 2 ). (2.7) 

3p 

Due to its computational complexity, the Boltzmann collision operator Qb(/) is often replaced in 
applications by simpler operators, like the BGK operator which substitutes the binary interactions 
with a relaxation towards the equilibrium of the form [3] 

Q B GK(f) = M (M[f] - f), (2.8) 

where p > 0. The validity of this operator in describing the physics of non equilibrium phenomena 
has been the subject of many papers in the past [9]. 

2.1. Fluid-limit and asymptotic-preserving methods. If we rescale the space and time 
variables in (12.11) as 



x' = ex, t' = et, (2.9) 
and omit the primes to keep notation simple, we obtain 

9tf + v>V*f=~Q(f) (2.10) 

where e is the Knudsen number a non dimensional quantity directly proportional to the mean free 
pat between particles. 
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Fig. 2.1. P £ is the original singular perturbation problem and P^ t its numerical approximation characterized 
by a discretization parameter At. The asymptotic-preserving (AP) property corresponds to the request that P^ t is 
a consistent discretization of P° as e — ¥ independently of At. 



Now integrating (|2.10l) against the collision invariants in the velocity space leads to the follow- 
ing set of non closed conservations laws 

d t (<f>f)+V x (v<t>f)=0. (2.11) 



Close to fluid regimes, the mean free path between two collisions is very small. In this situation, 
passing to the limit e — > we formally obtain Q(f, f) = from (|2.10p and so / = M[f\. Thus, at 
least formally, we recover the closed hyperbolic system of compressible Euler equations 

d t U + V x ■ F{U) = (2.12) 

with 

U = (4>M[f}) = (p,pu,E) T , 



F{U) = (v(f>M[f]) = (pu, qu®u+ pi, Eu + pu) T , p = pT, 

where / is the identity matrix. Note that the above conclusions are independent on the particular 
choice of Q(f, f) provided it satisfies (I2.4[) and admits Maxwellian of the form (|2.6[) as local 
equilibrium functions. 

For small but non zero values of the Knudsen number, the evolution equation for the moments 
can be derived by the so-called Chapman-Enskog expansion or Hilbert expansion [5]. These ap- 
proaches originate the compressible Navier-Stokes equations as a second order approximation with 
respect to e to the solution of the Boltzmann equation. In this case, however, the particular choice 
of the collision operator influences the structure of the limiting Navier-Stokes system. 

The construction of numerical schemes which are capable to capture the fluid-limit just de- 
scribed is closely connected with the notion of asymptotic-preserving schemes (see Figure 12. ip . 
Here, in agreement with [26J [30] we give the following definition of asymptotic preserving methods 
for equation (|2.1[) 

Definition 2.1. A consistent time discretization method for \2. 1\) of step size At is asymptotic 
preserving (AP) if, independently of the initial data and of the stepsize At, in the limit e — > 
becomes a consistent time discretization method for the reduced system I2.1ty) . 
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This definition does not imply that the scheme preserves the order of accuracy in time in the 
stiff limit s — > 0. In the latter case we will say that the scheme is asymptotically accurate (AA). 

In the sequel we will consider the development of asymptotic preserving and asymptotically 
accurate schemes using the general setting of IMEX Runge-Kutta methods. 

3. IMEX schemes for kinetic equations. First we introduce the general formulation of 
IMEX schemes for kinetic equations together with some preliminary definitions. 
An IMEX schemes applied to a kinetic equation of the form 

d t f + v-V x f = -Q(f) (3.1) 

£ 

reads [30] 

i— 1 v -. 

F® =f n - At^aijv ■ V x FW + At^Oij— Q(F^) (3.2) 

3=1 3=1 £ 

1 

f l+1 = f n -AtJ2 mv ■ V X F« + AtJ2 m-Q(F^ ). (3.3) 

i=l i=l 

The matrices A — (ay), a,-j = for j > i and A = (a y -) are v x v matrices such that the resulting 
scheme is explicit in v ■ V x /, and implicit in Q(f). In general, an IMEX Runge-Kutta scheme, 
is characterized by the above defined two matrices and the coefficient vectors w — (wi, .., w„) T , 
w = (u>i, .., w v ) T . Since computational efficiency in the case of kinetic equations is of paramount 
importance, in the sequel we restrict our analysis to diagonally implicit Runge-Kutta (DIRK) 
schemes for the source terms (ay = 0, for j > i). In fact, the use of a DIRK scheme is enough 
to assure that the transport term v ■ V x / is always evaluated explicitly. The type of schemes 
introduced can be represented with a compact notation by a double Butcher tableau 



c 


A 




~T 




W 



where the coefficients c and c are given by the usual relation 



c t ^^a lj , c i = ^a ij . (3.4) 

3=1 3=1 

Using vector notations the schemes can be written in compact form 

F = f n e + AtA L(F) + —AQ(F) (3.5) 

£ 

f n+1 = r + Atw T L(F) + —w T Q(F), (3.6) 

£ 

where e = (1, 1, .., 1) T el", F — {F^\ . . . , F^f , Q(F) = (Q(F^), Q(F^)) T and L(F) = 
(L{FW '),..., with L(F«) = -v-V x F^. 

We refer to [H [H [3D] for more details on the order conditions for IMEX schemes. Let 
us remark that IMEX schemes are a particular case of additive Runge-Kutta methods and so the 
order conditions can be derived as a generalization of the notion of Butcher tree [23j . In particular, 
under the assumptions c = c and w — w, mixed order conditions are automatically satisfied up to 
third order. 
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Before stating the main properties concerning asymptotic preservation it is useful to charac- 
terize the different IMEX schemes we will consider in the sequel accordingly to the structure of 
the DIRK method. Following [3] we have 

Definition 3.1. We call an IMEX-RK method of type A (see T3Uj ) if the matrix A £ R UXIJ 
is invertible, or equivalently ^ 0, i = 1, . . . , v. We call an IMEX-RK method of type CK (see 
J3j) if the matrix A can be written as 



a 



with a = (a.21, . . . , a v \) T £ R^^ 1 ) and the submatrix A £ M^ -1 ^ x ( v ~^ invertible, or equivalently 
an 7^ 0, i = 2, . . . , v. In the special case a — 0, wi = the scheme is said to be of type ARS (see 
jiy) and the DIRK method is reducible to a method using v — 1 stages. 

We will also make use of the following representation of the matrix A in the explicit Runge- 
Kutta method 



A= I ~ , (3.8) 




where 5 = (a 2 i, . . . , a uX ) T £ W^ 1 and A £ W~ lxi '~ 1 . 

The following definition will be also useful to characterize the properties of the methods in the 
sequel. 

Definition 3.2. We call an IMEX-RK method implicitly stiffly accurate (ISA) if the corre- 
sponding DIRK method is stiffly accurate, namely 

a vl = Wi, i ■=!,... ,v. (3.9) 

If in addition the explicit methods satisfies 

a^i^Wi, i = l,...,v (3.10) 

the IMEX-RK method is said to be globally stiffly accurate (GSA) or simply stiffly accurate. 
Some remarks are in order. 
Remark 1. 

• For type A IMEX schemes we have an ^ and an = thus c\ ^= c\ and we cannot 
assume the simplifying condition c = c. 

• For IMEX schemes of type A or type CK the GSA property implies w u — and w v ^= 
thus we cannot assume the simplifying condition w — w. Note that for GSA schemes the 
numerical solution is the same as the last stage value, namely f n+1 = F^ u \ 

• From the observations above it is clear that order conditions for GSA type A IMEX schemes 
are particularly restrictive since c ^ c and w w. 

3.1. Asymptotic preserving IMEX schemes. We give now conditions for an IMEX 
scheme to satisfy asymptotic preservation and asymptotic accuracy. Here we do not consider 
the computational challenges related to the inversion of the implicit collision operator Q(f). We 
will focus on these aspects in the second part of the paper. 

We can state the following theorem which show that type A IMEX schemes are asymptotic 
preserving and asymptotically accurate. 

Theorem 3.3. If the IMEX method is of type A then in the limit e -} 0, scheme i3.5\) - t3.6\) 
becomes the explicit RK scheme characterized by (A, w, c) applied to the limit Euler system 12.1 Sty . 
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Proof 

To prove the above theorem let us first multiply the IMEX method (|3.5p - (|3.6|) by the collision 
invariants 4>(v) = 1, v, v 2 and integrate the result in velocity space. We obtain the explicit Runge- 
Kutta methods applied to the moment system (|2.11[) 

(cfrF) = (0f e) + Ail(^(F)> (3.11) 
W n+1 ) = W n ) + Mw T (4,L(F)). (3.12) 



Now let us rewrite equation (|3.5[) in the following form 

eF = eFe + eAtAL(F) + AtAQ(F). (3.13) 
Since A is invertible we can solve for Q(F) to get 

AtQ(F) = eA- 1 (f - f n e - AtAL{F)^j . (3.14) 

As e — > we get 

AtQ(F) = F = M[F]. (3.15) 
Thus (j3.lip ~ p.12p becomes the explicit Runge-Kutta method applied to the limiting Euler system 

(232) 

U = U n e - AtA\7 x ■ F(U) (3.16) 

jjn+l _ jjn + Atli T Wx . (3 

where U = . . . , U^) T , F{U) = (.F(W (1) ), • ■ • ,F(U^)) T , U® = {(j>M[F^\) and F{U^) = 

((f>L(M[F^])). 

To conclude the proof we must be able to pass to the limit e — > in the numerical solution. 
This is in fact possible since for type A IMEX schemes the numerical solution is independent on 
e. In fact, we have 

= f l + Atw T L(F) + — — Q(F) (3.18) 

£ 

which thanks to (|3 . 14[) yields 

f n+1 = f n + Atw T L(F) + w T A- 1 (f - f n e - AtAL(F)^ 

and finally 

f n+l = f n ^ _ W T A -1^ + M _ ^-l^j L(F) + / ^-1 R (3^9) 



□ 

An important property of the schemes is that in the limit e — > the distribution function 
is projected over the equilibrium f n+1 — > M[f n+1 ]. From (I3.19P it is clear that this property is 
achieved if the following conditions are satisfied 



w T A~ 1 e=l, w T ^w T A- 1 A, w T A~ l M[F] = M[f n+l ]. (3.20) 
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The first condition corresponds to the classical L-stability requirement of the DIRK method [23 . 
Since the third condition depends on the stage values vector F the only possibility to satisfy p. 201) 
is that the IMEX scheme is GSA. In this case, in fact, we have 

W T A~ 1 = (0, . . . , 0, 1) T , M[F V ] = M[f n+l ), 

Thus, we can state 

Theorem 3.4. If the IMEX scheme is of type A and GSA then 

lim/" +1 = M[/" +1 l. (3.21) 

As observed, the request that the matrix A is invertible can be quite restrictive for high order 
methods since we cannot assume the simplifying condition c = c. However, under additional 
hypothesis, we can obtain schemes which are asymptotic preserving and asymptotically accurate 
even when the first row of A contains only zeros. 

In order to do this we first introduce the notion of initial data consistent with the limit problem. 

Definition 3.5. The initial data for equation h3.1\) are said consistent or well prepared if 



f (x,v)=M{f (x,v)]+g s (x,v), lim g s (x, v) = 0. (3.22) 

We can now state the following 

Theorem 3.6. // the IMEX scheme is of type CK and GSA then for consistent initial data 
in the limit e — > scheme Il3.5\) - [3~b]) becomes the explicit RK scheme characterized by (A,w,c) 
applied to the limit Euler system i2.12\) . 

Proof 

First we rewrite F = (^F^, F^j , w — (wi,w) T , w = (wi 1 w) T with F,w,w € IR^ -1 and the 

matrix A of the explicit tableau as in (|3.8[) . 

Then for type CK the IMEX schemes (|3.5|) - p.6j) can be written as 



F = f n e + AtiL(F^) + AtAL(F) + ^aQ(F^) + -^AQ(F), 



(3.23) 



fn +i = p + AMuZ^ 1 )) + Atw T L(F) + wi—Q(fW) + w T — Q{F) (3.24) 

£ £ 



where e = (1, .., 1) € W . The corresponding moment scheme yields 

= (4>f n ) 

(4>F) = (4>f n e) + At~a(4>L{F^)) + AtA{(f>L(F)) 
W n+1 ) = {4>D + Am(cf>L(F^)) + Atti T {4>L{F)). (3.26) 



(3.25) 



Solving now, the second equation in (I3.23[) for Q(F) we get 

AtQ(F) = eA- 1 \P - f n e - AtAL(F) - AtaL(f n )] - AtA^aQif 71 ). (3.27) 
As £ -> we obtain 

Q(F) = -A- l aQ{r). (3.28) 

Since /" = M[f n ] + g e in the limit e -> we have /" = M[f% Q(f") = 0, and ([3~28l) reduces 
to Q(F) = which implies F — M[F}. Moreover, since the scheme is GSA, we also have f n+1 = 
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M[f n+1 ] and thus at the next time step the initial value remains consistent and the moments 
system (|3.25[) - (|3.26p corresponds to the explicit Runge-Kutta methods for the Euler equations. 
□ 

If we now remove the assumptions of GSA and consistent initial data, the numerical solution 
takes the form 

/' 1+1 = /" (l - w T A- 1 ej + At (wi - tFi^fi) L(f n ) + At (ii T - w T A- 1 Tj L(F) 
+ w T A- 1 F+^(w 1 ~w T A- 1 a)Q(r)- 



Thus in order to be able to pass to the limit e — > we require the additional condition 

wi = w T A- 1 a. (3.30) 

The above equation is obviously satisfied by IMEX schemes of type ARS since w\ = and a = 0. 
In addition, if the IMEX scheme is of type ARS equation (|3.28j) reduces to Q(F) = which implies 
F = M[F]. Thus, except for the initial layer effect given by the O(At) term L(f n ) in (pT2"5)) - 
P-26P IMEX scheme of type ARS, even if not necessarily asymptotically accurate, are asymptotic 
preserving without further assumptions. The additional requirement wi = in this case suffices 
to give an asymptotically accurate method. 

Finally, we can state the following results on the asymptotic behavior of the numerical solution 

Theorem 3.7. If the IMEX scheme is of type CK and GSA then 

lim/" +1 = M[f n+1 ], (3.31) 

if one of the following conditions is satisfied 
(a) the initial data is consistent; 
(b) 

elA- l a = 0, (3.32) 

where e v = (0, . . . , 0, 1) T € M'^ 1 . 
Note that (|3.32|) is true for ARS type schemes. We end this paragraph with some remarks 
Remark 2. 

• In case the IMEX methods are not GSA the numerical solution may originate a final layer 
effect which leads to reduction of accuracy in the kinetic variable / n+ . Similarly for type 
CK IMEX schemes, in case the initial data is not consistent, the solution may exhibit an 
initial layer which gives rise to a reduction of accuracy in the numerical solution. These 
phenomena can be cured using smaller time steps or extrapolation techniques only for the 
very last, respectively first, step of the computation. 

• To avoid the loss of accuracy in the kinetic variable of non GSA IMEX methods it is also 
possible to impose the additional conditions up to third order J^j 

b T A- 1 c=l, b T A- 1 <? = l, b T 'A- 1 Ac= 1/2, (3.33) 

where c = Ae. For the type CK (consequently for the type ARS) we can rewrite these 
braic order conditions replacing A -1 with A . 



3.2. Properties of IMEX schemes for relaxation operators. In this paragraph we 
consider the particular case of BGK relaxation operators of the form QsGif(/) = f-i^lf) — f). 
A fundamental property of the IMEX scheme (|3.5|) - ()3.6j) applied to relaxation operators is that it 
can be solved explicitly despite the nonlinearity of M[f]. 
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To understand this let us remark that since the implicit scheme is a DIRK method equation 
(|3.2[) takes the form 



-l 

F (i) =f n + AtY,a ij L(F^)+AtJ2a ij ^(M[F^}-F (j '>) (3.34) 



where the only implicit term is the diagonal factor M[F^} — F^' in which M[F^] depends only 
on the moments (<f>F^ ) . If we now integrate equation (I3.34[) against the collision invariants thanks 
to the conservations (|2.4[) we obtain 

i-l 

(0F«) = (</>/") + At^a^L(F«)) (3.35) 
i=i 

which corresponds to the stages of the explicit Runge-Kutta method applied to the moment system 
(|2~TT|) . Thus (0F ( ^), and so MtFW], can be explicitly evaluated and system (|3T34|) is explicitly 
solvable. This property has been used, for example, in [3T| to implement efficiently IMEX methods 
for the BGK system. 

Monotonicity properties for Runge-Kutta methods have been studied by several authors in the 
recent past Q21 [22] ■ The same properties have been analyzed in the specific case of additive Runge 
Kutta schemes in [24]. Here, we restrict our study to the space homogeneous case (£(/) = 0) leaving 
to a future research the analysis in the non homogeneous case. Even if this case is particularly 
simple, we report briefly some results, since they are the basis of the analysis for the Boltzmann 
equation that we will discuss in the next section. Note, however, as shown in the above cited 
papers, that monotonicity conditions in the non homogeneous case are rather restrictive. Thus, 
if positivity is strictly demanded, as for instance in the case of Monte Carlo methods, splitting 
strategies are normally preferable |16j . 

In the homogeneous case the method reduces to a simple application of the implicit scheme 
and reads 

F = f n e+^A(M[f n ]e- F), (3.36) 
p+i = fn + w T(^i( M ^ e _ F ^ (3 37) 



where now the local equilibrium M[F) — M[f n ]e is independent of time since (4>f n+1 ) - 
Let us define z — fiAt/e and solve for F, we get 

F = (I + zA)- 1 (f n e + zAM[f n ]e) (3.38) 

f n+1 = (1 - w T A- 1 e)f l + w T A- 1 F 1 (3.39) 



where we used the fact that 



,(M[f n ]e-F) = A- 1 [F-f n e}. (3.40) 



Thus we have the following 

Proposition 3.8. Sufficient conditions to guarantee that f n+1 > when f n > in K3.36]) - 
( pO?P are that 

{I + zAy X e > 0, (I + zAy 1 Ae > 0, (3.41) 

l-w T A- 1 e>0, w T A^ 1 > 0. (3.42) 
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Note that conditions (|3.41l) - (|3.42l) must be interpreted component by component. In particular, 
(|3.4ip depend on z and (|3.42p require A to be invertible. For ISA schemes conditions reduce to 
(|3.4ip . The maximum value of z > such that conditions (13.411) are satisfied for z g [0,z] is 
usually referred to as the radius of absolute monotonicity of the scheme [18] . 
Let us remark that since 

(F - M[P]e) = (I + zA)- l {f n ~ M[f n ])e. (3.43) 

the numerical solution can be written as 

f n+l = R{z).f n + (1 - R(z))M[f n ], (3.44) 

where R(z) = 1 — zw T (I + zA)~ 1 e is the stability function of the DIRK method [23]. Note that 
(|3.44|) has the same structure of the exact solution to our problem where R(z) is an approximation 
of exp(— z). Then, we can state 

PROPOSITION 3.9. A sufficient condition to guarantee that f n+1 > when /" > in \3.36}) - 
|Q7| ) is that < R(z) < 1. 

Clearly the above proposition defines a subset of the absolute stability region of the method [23"] 
characterized by |-R(z)| < 1- Moreover, being based on a convexity argument, it gives a sufficient 
condition for the entropic property of the scheme, namely H(f n+1 ) < H(f n ) (see (|2.5[) V 

Remark 3. In this simple situation the L-stability of the implicit methods implies the AP 
property. In fact if R(z) -> R(oo) as z -> oo, then R(oo) = yields f n+1 — ¥ M[f n ] as s — ¥ 0. A 
recursive computation also shows that 

f n+1 = R(z) n+1 f + (1 - i?(z) n+1 )M[/ ], (3.45) 

since M[f n ] = AI[fo], Vrt. Therefore f n+1 — > M[fo] as n — > oo provided \R{z)\ < 1. Hence if 
< |-R(oo) | < 1 we have the weaker AP property f n+1 — > M[fo] as e — > 0, n — > oo. 

4. Penalized IMEX schemes for the Boltzmann equation. In the sequel, we will de- 
scribe how to modify the IMEX Runge-Kutta strategy in the challenging case of the Boltzmann 
equation (|2.1|) where the collisions are described by the operator Q(f) = Qb(/) introduced in 

Although the approach just described remains formally valid, in practice a major difficulty 
concerns the need to solve the system of nonlinear equations originated by the application of the 
DIRK method to the collision operator Qb(/)- The computational cost of such integral operator, 
characterized by a five fold nonlinear integral which depends on the seven dimensional space 
(x,v,t), is extremely high and makes it prohibitive the use of iterative solvers. Thus, to overcome 
this difficulty, the idea is to reformulate the collision part using a suitable penalization term. This 
idea, using a BGK model as penalization term, has been introduced recently in [TH]. We refer 
to [16j for an extension of this approach to exponential Runge-Kutta techniques. 

4.1. Penalization of the collision integral. The key idea of penalization techniques for 
the Boltzmann collision operator in stiff regimes is to observe that in such regimes the solution / 
is close to the Maxwellian equilibrium M[f]. Due to its computational and analytical complexity a 
large variety of simplified approximations to the Boltzmann collision operator which are valid close 
to equilibrium have been derived in the literature [3, 9, 25 . We mention here the BGK approxima- 
tion [3] QbgkU) — f*{M[f] — /) already discussed, which is capable to capture the leading order 
term, i.e. the compressible Euler equations. We mention also the ES-BGK relaxation approxi- 
mation [25] Q E s(f) = K M e[f\ ~ /) where a modified equilibrium M s [f], M [f] = M[f] is used, 
such that the equation has the advantage of matching the 0(e) expansion, i.e. the compressible 
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Navier-Stokes system. Finally, another possibility is to consider the linear Boltzmann equation in 
the form [5] 

QlbU) = I B{\v-vA,n){f(v')M[f]{v'*)-f{v)M[f}{v*)}dv*dn. (4.1) 

A fundamental property shared by all the previous operators is that their kernel is spanned by 
the local Maxwellian equilibrium M[f], so that the e — > asymptotic behavior of the original 
Boltzmann equation is preserved. 

Let us now denote with Qp(f) a general operator which will be used to penalize the original 
Boltzmann operator Qs(/). As discussed, the characteristics of Qp{f) are to be computable 
and invertible at a low computational cost and that it preserves the local equilibrium, namely 
Q P (f) = implies / = M[f]. 

We will then rewrite the collision operator in the form 

QbU) = (Qs(f) - Qp(f)) + Qp(f) = Gp(f) + Qp(f), (4.2) 
where by construction (cf)Gp(f)) = 0, and the corresponding kinetic equation reads 

dtf + v-W x f = -Gp(f) + -Q P (f). (4.3) 

e e 

Now, the idea is to use a numerical scheme in which only the simpler operator Qp(f) is treated 
implicitly, while the term Gp{f) describing the deviations of the true Boltzmann operator Qb(J) 
from the simplified operator Qp(f) and the convection term V x / are treated explicitly. This 
however, as we will see, introduces some additional stability requirements in order for the IMEX 
schemes to preserve the asymptotic behavior of the equation. 

Let us remark that the above penalization method in the case Qp(f) = Qp,GK(f) admits an 
interesting interpretation that relates the previous approach with decomposition methods |151 1 141 
[27] . First, we observe that we can rewrite Qb(J) as 

Qb(J) = P(f) - tf, (4.4) 

where P(f) = Qb(I) + an d /i > is a constant such that P(f) > 0. Typically [i is an estimate 
of the largest rate of the negative term in the Boltzmann operator 

M > / B(\v — v*\, n)/(u«) du* dn. (4.5) 

JR 3 xS 2 

By construction, the following property is verified by the operator P(f, f) 

-^P{fJ)) = W) = U. (4.6) 

Thus P(f, f)/fi> is a density function and we can consider the following decomposition 

P(fJ)/H = M[f]+g, (4.7) 

where the function g = P(f,f)/fi — M[f] represents the non equilibrium part of the distribution 
function and from the definition above it follows that g is in general non positive. Moreover since 
P{f,f)/n and M[f] share the first three moments we have ((fig) = 0. 
Thus the collision operator can be rewritten in the form 

Q B (f, f) = nM[f] +ng- nf = ng + fi(M[f] - f) = G BGK (f) + QbgkU), (4.8) 

which is the same as (|4.2[) using the BGK model as penalizer. As we will see, thanks to the 
derivation above, the use of the BGK model is particularly interesting to study the monotonicity 
properties of the penalized IMEX schemes. 
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4.2. Asymptotic preserving penalized IMEX schemes. We can now introduce the gen- 
eral class of penalized IMEX Runge-Kutta schemes for the Boltzmann equation in the form 

F« = /™ + A^a y (-G P (F^) - v -V X F(A + Atf2<*ij-Qp{F lj) ) ( 4 -9) 

i=l V£ ' 3=1 £ 

f" +1 = f n + AtJ^m (±G p (fU)-vV x f(A +AtY l w i jQp(F(% (4.10) 

i— 1 ^ ' 

or equivalently, using vector notation as 



F = f n e + AtA ( -Gp(F) + L(F) ) + AtA-Q P (F) (4.11) 



fn+ i = f n + At{b i \-G P (F) + L(F)\ + Atw 1 -Qp(F), (4.12) 

where G P (F) = (G P (F^), G P {F^)) T . 

We want now to derive conditions for the penalized IMEX schemes to be AP and asymptotically 
accurate. As we will see, even for type A IMEX schemes, additional conditions are required to 
achieve asymptotic preservation. 

Theorem 4.1. If the penalized IMEX method is of type A and satisfies 

w T = w T A~ 1 A, (4.13) 

then in the limit e — > 0, scheme ^.11^ - ^.12 ) becomes the explicit RK scheme characterized by 
(A,w,c) applied to the limit Euler system \2.12\l . 
Proof 

First let us observe that if we multiply the penalized IMEX method (|4.11l) - (|4.12l) by the 
collision invariants <j)(v) — l,v,v 2 and integrate the result in velocity space we obtain again the 
explicit Runge-Kutta methods applied to the moment system (|2.11l) 

(cf>F) = {4>f n e) + AtA{cf>L(F)) (4.14) 
W n+1 ) - (4>f n ) + Atw T (<t>L(F)). (4.15) 

Now let us rewrite equation (|3.5|) in the following form 

eF = eFe + eAtA \-G P {F) + L{F) ] + AtAQ P (F). (4.16) 



Since A is invertible we can solve for Qp(F) to get 

AtQp(F) = eA- 1 (f - f n e - AtA (^G P (F) + L(F)X\ . (4.17) 

As e — > we get 

Q P {F) = -A- X AG P (F). (4.18) 

Since A~ 1 A is lower triangular with diagonal elements equal to zero we get 

Q P (F (i) )=0 F®=M[F®], i = l,...,u. (4.19) 
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Thus (|4.14j) - (l4.15l) becomes the explicit Runge-Kutta method applied to the limiting Euler system 

(E32). 

Unfortunately, now the numerical solution still depends on e. In fact, inserting (|4.17l) into 
(|4.12p we obtain 

jn+i = f n + am t (~G P {F) + L{F) \ + w T A- 1 (f - f n e - AtA (^G P (F) + L(F) 
which gives 

f n+i = fl ^ _ w T A -i e j + At _ ^ A -i^j (±G P (F) + L(F)) + w T A' 1 F. 

Using now the assumption (|4.13|) we get 

p+i = p ^ _ w T A -i e ^ + w t a ~i f ^ ( 4 20) 

which permits to pass to the limit e — > in (|4.11[) - (|4.12[) . 

□ 

Note that condition (|4.13|) is automatically satisfied if the IMEX scheme is GSA. In this case 
it is easy to show that we have 

Theorem 4.2. // the penalized IMEX scheme is of type A and GSA then 

lim/" +1 = M[f n+1 }. (4.21) 

We consider now the case of penalized IMEX schemes of type CK. We can state an analogous 
result of Theorem 13.61 for standard IMEX schemes of type CK. Some relevant differences will be 
pointed out at the end of the proof. We have the following 

Theorem 4.3. // the penalized IMEX scheme is of type CK and GSA then for consistent 
initial data in the limit e — > scheme \^.11^ - {4.12^ becomes the explicit RK scheme characterized 
by (A,ib,c) applied to the limit Euler system \2.12\l . 

Proof 

Using the same notations as in Theorem [376] the penalized IMEX schemes (|4.11ll - (|4.12p can be 
written as 

pi 1 ) — fn 

F = f n e + Ata (-G P (F^) + L(F^)) + AtA (-G P (F) + L{F)) (4.22) 



e I \e 



+ ^aQ P (F^) + —AQ P (F), 

e e 

f n+1 =f n + Atwi (-G P (F^) + L(fW)) + Atw T (-G P {F) + L{F) 



e I \e 



■wr—Q P {F^)+w T —Q P {F), 

e e 



(4.23) 



whereas the corresponding moment scheme is unchanged and reads 

(d>F^) = {4>n 

= (cf>f n e) + Ata{(/>L(F^)) + Atl{cf>L{F)) 
= (</>/") + Atw^LiF^)) + Atti T (0L(F)). (4.25) 
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(4.24) 



F - f n e - AtA ( -Gp(F) + L(F) I - Ata I -G P (f n ) + L(f n ) 



Solving now, the second equation in (|4.22[) for Qp(F) we get 
AtQp(F) = eA- 1 

' ' (4-26) 
- AtA-'aQ P {f n ). 

As £ — > we obtain 

Qp(F) = -i- 1 (IGp(F) + ~aG P {P)) - A^aQpif 1 ). (4.27) 

Now since /" is consistent in the limit e^Owe have /" = M[f n ], Qp(f") = 0, G P (f") = and 
(R~2"T)) reduces to 

Qp{F) = -A- l AGp{F) (4.28) 

which thanks to the fact that A -1 A is lower triangular with zero diagonal elements implies F — 
M[F]. Moreover, since the scheme is GSA, we also have f n+1 = M[f n+1 ] and thus at the next 
time step the initial value remains consistent and the moments system (I4.24p ~ d4.25p corresponds 
to the explicit Runge-Kutta methods for the Euler equations. 
□ 

Note that even for ARS type IMEX schemes with a = the diagonal elements of the right 
hand side of equation (|4.27p do not vanish and thus asymptotic preservation is not achieved unless 
the initial data is consistent or the following additional condition holds true 

~A- 1 i = 0. (4.29) 

In this latter case, ARS type schemes are asymptotically accurate if w\ = 0. 

Let us consider now the behavior of the penalized IMEX schemes of type CK in case the GSA 
and the consistent initial data assumptions are removed. The numerical solution takes the form 

r +i - r (\ - uFA- 1 ^ + At (w x - w T A- iz a) Qc P (r ) + + ^a^p 

+ At - uFA- 1 !) Qgp(F) + L(F) \ +^(wx- ^A^aj Q P {f n ), 
and thus we need the further assumptions 

W! = w T A~% w T = w T A- l A, wi = w T A~ l a, (4.30) 

in order to pass to the limit e — > 0. Clearly the assumptions above are satisfied by GSA schemes. 

We end the paragraph by emphasizing the following result on the asymptotic behavior of the 
numerical solution 

Theorem 4.4. If the penalized IMEX scheme is of type CK and GSA then 

lim/" +1 = M[f n+1 ], (4.31) 

if one of the following conditions is satisfied 
(a) the initial data is consistent; 
(b) 

elA- 1 A = 0, elA- 1 a = 0, e^A^a = 0, (4.32) 
where e v = (0,...,0, 1) T € R" -1 ; 

(c) 

A- 1 a = 0, A- 1 a = 0. (4.33) 
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4.3. Properties of penalized IMEX-BGK schemes for the Boltzmann equation. In 

this last paragraph, we restrict to the particular case where the penalization term is given by 
a BGK relaxation operator. Exactly as for standard IMEX schemes a fundamental property of 
equations (I4.1ip - (I4.12[) with Qp(f) = fi(M[f] — /) is that they can be solved explicitly. 
In fact, since the implicit scheme is a DIRK method the stage values take the form 

i-1 / -, \ i 

F (i) =f n + AtJ2aij [L(F^) + -G P (F j )) +AtJ2^-(M[F^}-F^) (4.34) 
3=1 ^ 6 ' j=i £ 

where the only implicit term is the diagonal factor M[F^'] — F^' in which M[F^] depends only on 
the moments (<f>F^). If we now integrate equation (|4.34|) agains the collision invariants thanks to 
the conservations (|2.4p we obtain again the moment scheme (|3.35[) . Thus (<f>F^), and so M[F^], 
can be explicitly evaluated and system (|4.34[) is explicitly solvable. This property was at the basis 
of the penalized IMEX schemes recently proposed in [19] . 

Next we focus on the monotonicity properties of the penalized IMEX schemes. It appears 
clear that in space non homogeneous situations due to the penalization approach it is extremely 
difficult to achieve positivity of the numerical solution. For this reason we restrict to the space 
homogeneous case, which in this case however involves the full penalized IMEX method. As usual, 
in combination with time splitting strategies, these results permit to recover positivity even in the 
non homogeneous case if required. This is of paramount importance, for instance, in the case of 
Monte Carlo methods [29] , 

In the homogeneous case the method takes the form 

F = f n e + t^-A ( - M[f n ]e) + --^-A {M[f n ]e - F) , (4.35) 
e \ n ye 

f l+1 = r + w T — ( — - M[f n ]e) + (M[r] e - F) , (4.36) 



where again the local equilibrium M[F] — M[f n ]e is independent of time since (4>f n+1 ) = (4>f n )- 
In (14351) - KM we used the fact that G P (F) = P(F) - nM[f n }e with P(F) = Q B {F) + ^F > 
and {<f>P(F)) = (tpnF). 

Let us define z — fiAt/e and solve for F, we get 

F=(I + zA)- 1 (f n e + z(A - A)M[f n ]e + zA^---) (4.37) 



= (1 - w T A- 1 e)f n + w T A- 1 F + z(w T - w T A' 1 A) { ^ - M[f n ]e ) , 1 4.:-!<s , 



where we used the identity 

z(M[f n ]e — F) = A- 1 



F- /"< - .-..1 ( ^1-M[f n \e 



(4.39) 



We can thus state the following 

PROPOSITION 4.5. Sufficient conditions to guarantee that f n+1 > when f n > in U-S5\ )- 
£p]6\ ! are that 

(I + zAy 1 e>0, (I + zAy 1 (A-A)e > 0, (I + zA)' 1 A>0 (4.40) 
l-w T A- 1 e>0, w T ^ _1 >0, w T ^w T A- 1 A. (4.41) 
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Conditions (I4.40j) - (14.41 j) applies to a general penalized IMEX method and must be interpreted 
component by component. Clearly for GSA schemes conditions reduce to (I4.40[) . Since (14.401) - 
(14.41)) are originated from the penalizing strategy in the nonlinear case they differ from the classical 
conditions for absolute monotonicity of additive Runge-Kutta methods [24] . The range of values of 
z for which (I4.40|) - (I4.41|) hold true dehnes the absolute monotonicity region of the penalized IMEX 
method. Note that here we are interested in the behavior of the schemes outside the stability 
region of explicit methods, thus monotonicity is important for z > 1. 

Finally, in order to characterize the stability property of the penalized IMEX schemes we 
consider the linearized problem 

dtf=-G L (f) + ^(M[f]-f), (4.42) 
e e 

where Gl(/) = (A — /i)(M[/] — /). Here fi is an estimate of the unknown relaxation rate A > 0. 
We will further denote bya = A//i>0 the penalization parameter of the method. Clearly a = 1 
corresponds to the optimal choice for which the whole analysis reduces to the one performed in 
Section [3T2l for the BGK case. 

The penalized IMEX schemes take the form 

F = f n e + zA(a — 1) (F — M[f n ]e) + zA (M[f n ]e - F) , (4.43) 

f n+1 = P + w T z(a — 1)(F — M[f n ]e) + w T z (M[f n ]e - F) . (4.44) 

Now since 

(F - M[P]e) = (I+ z(A - (a - lJA))"^/" - M[/"])e. (4.45) 
the numerical solution becomes 

f n+1 = R(a, z)P + (1 - R(a, z))M[p], (4.46) 

where 

R{a, z) = l- z{w T - (a - l)w T )(/ + z{A - {a - l)i)) _1 e (4.47) 

plays the role of the stability function of the penalized IMEX method. In this case R(a, z) is 
an approximation of exp(— az) and we characterize the absolute stability region by the condition 
\R(a, z)\ < 1. Note that R(l, z) corresponds to the stability function of the implicit method. 
In particular, similar to the BGK case (|3.45l) . we recursively obtain 

P+ 1 = R(a, z) n+1 p + (1 - R(a, z)" +1 )M[/ ]. (4.48) 

This permits to consider a weaker notion of AP property when \R(a, oo)| < 1 where 

R(a, oo) = lim R(a, z). (4.49) 

In this case intact f n+1 — » M[fo) as n — > oo at a rate characterized by the damping constant 
|i?(oo,a)| < 1. Of course for an L-stable implicit method we have R(l, oo) = 0. From (|4.49[) an 
AP method corresponds to the requirement R(a, oo) — 0, Va > 0. A condition which essentially 
generalizes the notion of L-stability to penalized IMEX schemes. 

The sensitivity of i?(a, oo) to values of the penalization parameter a / 1 gives a measure of 
the robustness of the method with respect to the estimate fi of A in terms of capturing the correct 
asymptotic behavior. Thus we introduce the following 
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Table 5.1 

Summary of the properties of penalized IMEX GSA schemes 



Name 


AA 


AA-c 


AM 


R(a, oo) 


\R{a,oo)\ < 1 


ARS(1,1, 1) 


no 


yes 


z e [0, oo) 


a - 1 


aG (0,2) 


DP-ARS*(1,2,1) 


yes 


yes 


2 = 





a 6 (0, oo) 


DP-A (1,2, 1) 


yes 


yes 


z G [0, oo) 





a £ (0, oo) 


AKo(z, I, 1) 


no 


yes 


z = U 


(l — a)(a+4~f — 2a-f — 2~f 2 ) 


a G (U.o 1 4, 1.11 ( J 


DP-ARS(2,2,2) 


no 


yes 


z e [o, oo) 


(l-a)(a+47-2d7 — 27 2 ) 
2 7 2 


a G (0,2.288) 


JF-CK(2,3,2) 


no 


yes 


ze[o,2] 


2a 2 - 4a + 1 


a G (0,2) 


DP1-A*(2,4, 2) 


yes 


yes 


2 = 





a G (0, oo) 


DP2-A 2 (2,4,2) 


yes 


yes 


2 G [1, oo) 





a G (0, oo) 


ARS(4,4, 3) 


no 


yes 


2 = 


(l-a)(7a 3 -78a 2 +138a-38) 
18 


a G (0.13475,2) 


BPR-CK(3,5,3) 


no 


yes 


2 = 


4(2a 3 -5a 2 +2a) . , 

3 + 1 


ae(0,i^)U(i,i^)U(3,2) 



Definition 4.6. The values of a > for which \R(a, oo)| < 1 characterizes the weak AP 
range o/ i/ie penalized IMEX method. 

Remark 4. A^oie t/iat, achieve better penalization properties, the value of fj, can be taken 
space and time dependent in ^4-'M^ t° minimize the value of Gp(F^). Typically, for a given f, 
one should choose /i > such that it locally minimizes \\P(f) — [J.M[f]\\ in a suitable norm || • ||. 
Here we do not explore further this aspect. However, the general asymptotic preserving properties 
and considerations we have done remain valid also in this case. 

5. Numerical considerations and tests. We report in Table [Ol a summary of the proper- 
ties of some penalized IMEX schemes satisfying the GSA property developed in the literature. We 
use the notation NAME(^e, vi,p) where ve, vi are, respectively, the number of function evaluations 
of the explicit and the implicit methods and p is the combined order of the IMEX scheme. The field 
NAME of the schemes is composed by the initials of the authors and the scheme type. We con- 
sider the following methods: ARS(1, 1, 1) corresponding to forward-backward Euler, ARS(2, 2, 2) 
and ARS(4,4,3) from Sections 2.6 and 2.8 in [I], DP-ARS(2, 2, 2) same as ARS(2,2,2) but with 
7 = (2 + V2)/2, JF-CK(2,3,2) from (2.8) Section 2 in [19], BPR-CK(3, 5, 3) from the Appendix 
in g]. Schemes DP1-A(1, 2, 1), DP2-A 2 (1,2,1), DP-ARS(1, 2, 1) and DP-A(2,4,2) are reported in 
a separate Appendix at the end of the manuscript. With the acronym AA we denote the asymp- 
totic accurate property, with AA-c the same property but for a consistent initial data. We also 
report the region of absolute monotonicity (AM) in the nonlinear case accordingly to definition 
(|4.40[) - (14.41[) . the asymptotic behavior of the stability function (|4.49|) and the weak AP range of 
the schemes obtained from Definition 14.61 For schemes marked with an '*' the DIRK part has 
order p + 1. 

We emphasize that the computational cost of penalized IMEX schemes is characterized by 
the number of stages of the explicit method since, by construction, the implicit part is applied 
to the easy invertible term used for penalization. This is an important aspect for the practical 
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construction of the methods. 

Next we numerically measure the accuracy of second and third order methods in solving the 
full Boltzmann equation (|2.1[) - (|2.2I) on a periodic smooth solution. We emphasize that to our 
knowledge these are the first full third order implementations for the Boltzmann equation. 

The computation is performed on (x,v) € [0,1] x [— u max , v max ] 2 , with u max = 8. We use a 
3rd order WENO scheme for the space discretization [32] and a fast spectral method for solving 
the collision integral [28]. We take N v = 32 grid points in each velocity direction which provides 
enough accuracy in v to measure the convergence rates in time. The time step is fixed equal to 
At = Ax/(2v mayi ). The initial data is 

2 + sin(27re) cos(27rx) 3 + cos(27ra;) 
Qo(x) = , Uo{x) = , Tq{x) = . (5.1) 

We report the L\ norm of the error for the density for different values of the Knudsen number, i.e. 
e = 1CP 1 , e = 10~ 3 and e = 10~ 6 . We consider both equilibrium initial data fo(x,v) = M[fo] and 
non equilibrium initial data 

fo(x,v) = ™) ' - exp 2T oM +exp 2T oM . (5.2) 



(27rTo(a:)) 1 /2 2 

To measure the convergence rate we repeat the computation for an increasing number of grid points 
in space, N x = 128, 256, 512, 1024. 

Figure |5~T1 shows the results for the second and third order penalized IMEX schemes. Here in 
scheme DP2-A(2,4,2) we choose to maximize accuracy taking 7=1/3 outside the monotonicity 
region. As expected, all the schemes exhibit the prescribed order of convergence for equilibrium 
initial data while degradation of accuracy to first order is observed for schemes not asymptotically 
accurate. Schemes DP1-A(2, 4, 2) and DP2-A(2, 4, 2) preserves second order accuracy uniformly in 
all regimes independently of the initial data. Finally, we report in figure [5T21 the comparisons be- 
tween the convergence rate of the second order monotone schemes where the diagonal value is taken 
inside and outside the monotonicity region. Note how the large diagonal entry needed for mono- 
tonicity reduces the accuracy of the schemes for non equilibrium initial data in the intermediate 
regimes. 

6. Conclusions and further developments. We have studied the extension of implicit- 
explicit (IMEX) Runge-Kutta schemes to the numerical solution of nonlinear kinetic equations in 
stiff regimes. We derived sufficient conditions for an IMEX scheme to be asymptotic preserving and 
asymptotically accurate. These notions are of paramount importance in the practical computation 
of kinetic equations where regions with a large variation of the computational space and time scales 
are present. 

In the first part we focused on standard IMEX methods and studied their asymptotic behavior 
and monotonicity property. In the second part, following the approach introduced in [19] , we tackle 
the challenging case of the efficient computation of the Boltzmann collision term in stiff regimes. A 
penalization technique permits to derive new schemes, here referred to as penalized IMEX schemes, 
which under additional assumptions keep the asymptotic properties of the standard IMEX methods 
by avoiding the costly inversion of the collision term. Monotonicity is studied for the case where 
the penalization term is represented by a relaxation operator. 

For the future extensive numerical testing and applications of the schemes to other collisional 
kinetic equations, like the Landau-Fokkcr-Planck equation, is scheduled. 
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FlG. 5.1. Li error for the density q for different second and third order IMEX schemes. Left column equilibrium 
initial data, right column non equilibrium initial data. Top e = 10" 1 , center e = 10~ 3 , bottom e = 10 -6 . 



Appendix. In the sequel we report the Butcher tableu of some penalized IMEX schemes 
satisfying the AA condition. The guidelines in the design of the methods are the GSA property 
and the minimization of the number of explicit function evaluations. We also consider schemes 
for which the explicit part is strongly stability preserving |22| . As we will see since the above 
requirements lead to schemes using several implicit levels, it is natural to satisfy higher order 
conditions for the DIRK part. 
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Fig. 5.2. Li error for the density q for second order monotone schemes for different values of the diagonal 
entries in the monotone and non monotone regions. Non equilibrium initial data with e = 10 — 3 . 

Penalized first order schemes satisfying the A A property are easy to construct, we give here 
two examples, one absolutely monotone of type A, referred to as DP-A(1, 2, 1) 
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where 7 > 1/2, and a non monotone variant of type ARS that we refer to as DP-ARS(1, 2, 1) 
which satisfies the asymptotic preservation conditions (|4.32[) 
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where 5 = 7/(1— 7). For both first order schemes the choice 7 = (1 ± y/2)/2 yield second order 
accuracy for the DIRK part. 

Note that even if two levels are required for the penalized schemes to be GSA and achieve AA, 
they use just one explicit function evaluation. 

For GSA type A schemes it is easy to show that we need at least four levels to satisfy the 
second order accurate requirement and the GSA property. On the other hand four levels do not 
suffice to have third order accuracy and GSA property As examples we report below scheme 
DP1-A(2, 4, 2), non monotone, which achieves third order accuracy on the DIRK part 
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and a second order scheme DP2-A(2, 4, 2) 
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which is monotone on the interval z > {2^ 2 — 47 + for 7 > (2 + \/2)/2. In particular for 7 = 2 
monotonicity is obtained in [l,oo). From the accuracy viewpoint the choice 7 = 1/3 originates a 
DIRK method which satisfies the third order condition w T Ac =1/6. We denote by DP2-Ai(2, 4, 2) 
the scheme with 7 = 1/3 and by DP2-A 2 (2, 4, 2) the scheme with 7 = 2. 
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